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ABSTRACT 


Modeling of the dynamic vibration modes of a flexible structure can be achieved either by 
using a generalized coordinate for each mode considered in the simulation, or by 
discretizing the structure into a sufficiently large number of segments to provide the 
necessary modal accuracy. The accuracy and stability considerations in choosing 
appropriate numerical integration algorithms are different, depending on which 
modeling approach is utilized. In the generalized coordinate approach the frequency and 
shape of each mode is assumed to be known. The integration method should provide an 
accurate match to the modal frequency and damping, and should also exhibit sinusoidal 
transfer function errors which are acceptably small, especially for frequencies in the 
vicinity of the modal resonance. Since only those inodes considered necessary for the 
required simulation fidelity are included as generalized coordinates, integrator stability 
for modes of higher frequency does not become in an issue. 


On the other hand, when the discretized structure approach is used, high frequency modes 
not of interest to the simulation will nevertheless be present. In this case it is important 
that the integration method not only provide satisfactory characteristic root and transfer 
function accuracy for the lower modes of interest, but also provide stable solutions wit 
satisfactory damping for the higher modes which are not of interest. 


ir. iftio napet asymptotic formulas for the characteristic root errors as well as transfer 
function gain "and phase errors are presented for a number of traditional integration 
methods and for several new integration methods. Normalized stability regions in the 
Xh plane are compared for the various methods. In particular, it is shown that a modified 
form of Euler integration with root matching is an especially efficient method for 
simulating lightly-damped structural modes. The method has been used successfully for 
structural bending modes in the real-time simulation of missiles. Performance of this 
algorithm is compared with other special algorithms. Including the state-transition 
method. A predictor-corrector version of the modified Euler algorithm permits it to be 
extended to the simulation of nonlinear models of the type likely to be obtained when 
using the discretized structure approach. 

Performance of the different integration methods is also compared for integration step 
sizes larger than those for which the asymptotic formulas are valid. It is concluded that 
many traditional integration methods, such as RD-4, are not competitive in the 
simulation of lightly damped structures. 
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A Performance Comparison of Integration Algorithms in Simulating Flexible 

Structures 

R. M. Howe 

The University of Michigan, Ann Arbor, Michigan 
and 

Applied Dynamics International, Ann Arbor, Michigan 
ABSTRACT 

In this paper a number of integration algorithms, including several new methods, are 
considered for the simulation of flexible structures. The effectiveness of the different 
algorithms is assessed by considering the characteristic root errors which they produce, the 
sinusoidal transfer function gain and phase errors, the stability regions, and the execution times. 
The suitability of the various algorithms for simulations with real-time inputs is also noted. 
When the structural modes in a simulation are represented by generalized (normal) coordinates, 
the selection criteria for integration methods are somewhat different than the criteria when the 
structure is discretized into a sufficiently large number of segments to provide the necessary 
modal accuracy. In this paper asymptotic formulas for the characteristic root errors as well as 
transfer function gain and phase errors are presented for a number of traditional integration 
methods and for several new integration methods. Normalized stability regions in the X.h plane 
are compared for the various methods. In particular, it is shown that a modified form of Euler 
integration with root matching is an especially efficient method for simulating structural modes. 
The method has been used successfully for structural bending modes in the real-time simulation 
of missiles. A predictor version of the modified Euler algorithm permits it to be extended to the 
simulation of nonlinear models of the type likely to be obtained when using the discretized 
structure approach. 

1. Introduction 

Modeling of the dynamic vibration modes of a flexible structure can be achieved either 
by using a generalized coordinate for each mode considered in the simulation, or by discretizing 
the structure into a sufficiently large number of segments to provide the necessary modal 
accuracy. In this latter case the mathematical model for a flexible structure with N degrees of 
freedom has the following general form: 

M{q)q + C(q,q) + K(q) = F(t) (1) 

where q is an ^/-component position state vector, M{q) is the mass matrix, C(q,q) is the coriolis 
and centrifugal acceleration vector, K(q) is the elastic and gravity force vector, and F(t) is the 
external force vector. When the vibration modes of the structure are represented by normal 
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(generalized) coordinates, a coordinate x representing the time-varying amplitude of a given 
mode with undamped natural frequency (On and damping ratio £ obeys the equation 

x + 2Cco n x + (olx = cotyt) ( 2 ) 

Here 0(r) is the generalized force associated with the coordinate x. When a number of modes 
are present, there will in general also be terms in Eq. (2) which couple the mode of amplitude x 
with other structural modes. 

The accuracy and stability considerations in choosing appropriate numerical integration 
algorithms for solving differential equations of the type shown in (1) or (2) will be different. In 
the generalized coordinate approach of Eq. (2) the frequency and shape of each mode is 
assumed to be known. The integration method should provide an accurate match to the modal 
frequency and damping, and should also exhibit sinusoidal transfer function errors which are 
acceptably small, especially for frequencies in the vicinity of the modal resonance. Since only 
those modes considered necessary for the required simulation fidelity are included as 
generalized coordinates, integrator stability for higher frequency modes which are not of interest 
does not become an issue. 

On the other hand, when the discretized structure approach represented by Eq. (1) is 
used, high frequency modes which are unimportant in the simulation will nevertheless be 
present. In this case it is important that the integration method not only provide satisfactory 
characteristic root and transfer function accuracy for the lower modes of interest, but also 
provide stable solutions with satisfactory damping for the higher modes which are not of 
interest. 

In this paper asymptotic formulas for the characteristic root errors as well as transfer 
function gain and phase errors are presented for a number of traditional integration methods and 
for several new integration methods. Normalized stability regions in the Xh plane are compared 
for the various methods, where A is an eigenvalue asociated with the linearized perturbation 
equations of the structure and h is the integration step size. . In particular, it is shown that a 
modified form of Euler integration with root matching is an especially efficient method for 
simulating lightly-damped structural modes. The method has been used successfully for 
structural bending modes in the real-time simulation of missiles. Predictor versions of the 
modified Euler algorithm permit it to be extended to the simulation of nonlinear models of the 
type likely to be obtained when structures are represented by means of discretization. The 
stability regions in the Xh plane for the modified Euler methods are especially well suited to the 
requirements when using the discretized structure approach. 

2. Dynamic Error Measures for Integration Algorithms 

In comparing different integration methods for the simulation of flexible structures it is 
important to utilize meaningful performance measures which permit general conclusions to be 
drawn regarding the expected dynamic errors associated with each method. Our dynamic error 
analysis will be based on linearized perturbation equations derived from the original nonlinear 
differential equations used to model the structure. Thus we will assume that the system 
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eigenvalues are known, as well as the transfer functions relating specific input-output pairs. We 
will further assume that the simulation uses a fixed integration step size h. This is necessary in 
the case of a real-time simulation. It is likely to be true over a large number of steps even when 
a variable-step integration method is used in simulating a flexible structure. For linearized 
equations and a fixed integration step size we can apply the method of z transforms to anayze 
the dynamic errors resulting from specific integration algorithms [1,2]. There are two error 
measures which quite useful in predicting overall dynamic accuracy in the simulation. The first 
is the fractional error in each characteristic root (eigenvalue) of the digital simulation, defined as 


Fractional error in characteristic root = e A = 


where X is the characteristic root of the continuous system being simulated and X* is the 
equivalent characteristic root for the digital simulation. For the case of complex roots (of which 
there will be many conjugate pairs in the simulation of a flexible structure) it is more appropriate 
to determine the fractional error, e w , in root frequency and the damping ratio error, e^ Thus we 

define 


og- Qfr 


e { ={*-C 


(4) 


Here (O and (O d represent the frequencies of the digital and continuous system roots, 
respectively, while £* and £ represent the damping ratios for the digital and continuous system 

roots, respectively. 

The second dynamic error measure of significance is the fractional error in digital system 
transfer function for sinusoidal inputs of frequency (t). For any input-output pair let H(s) be the 
transfer function of the continuous system and H*(z) be the z transform of the digital system 
that results when a particular integration algorithm is used. Then the fractional error in 
sinusoidal transfer function is given by [3] 


H*(c JbJh ) 

H(ja>) 


1 = 


e M + J e A 


(5) 


For simulations of any reasonable accuracy the magnitude of this fractional error will be small 
compared with unity, in which case it is easily shown that the real part, eM, is equal 
approximately to the fractional error in gain and the imaginary part, e a, is equal to the phase 
error of the transfer function [3]. 

For any numerical integration algorithm the integrator transfer function for sinusoidal 
inputs of frequency CD can be written approximately as [3] 


H*(c jah ) = 


ja>h^ \ + e / (/'mfc)*j 


coh « 1 


(6) 
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where h is the integration step size. Since 1 /(jcoh) is the ideal integrator transfer function, it is 
apparent that the term ejijoh^ represents the integrator error. For Adams-Bashforth predictor 
and Adams-Moulton two-pass predictor-corrector algorithms of order 2, 3, and 4, integration 
methods that are candidates for simulation of flexible structures, the error coefficient e, and 
algorithm order k are listed in Table 1 . 


Table 1. Integrator Transfer Function Error Parameters for AB Predictor and 

AM Predictor Corrector Algorithms 


Method 

Error Coefficient, ej 

Algorithm Order, k 

AB-2 

5 

TT 

2 

AB-3 

3 

S’ 

3 

AB-4 

251 

720 

4 

AM-2 

1 

TT 

2 

AM-3 

l 

3* 

3 

AM-4 

19 

72U 

4 


In terms of e\ and k the following formula for ex, the fractional error in characteristic root as 
defined earlier in Eq. (3), can be derived [3]: 

A* X 

- T = - e ; (Xh) k , \Xh\« 1 ( 7 ) 

It is apparent that e\ is directly proportional to the integrator error coefficient, ej. For complex 
characteristic roots equivalent asymptotic formulas for the root frequency and damping errors 

c'bS H ' aS defmed in Eq ' (4) ’ Can be derived [3] ‘ As in Eq - (7) * the erTOrs Proportional to 

For digital simulation of a first order system with transfer function H(s) = \/{s-X) the 
fractional error in the transfer function for sinusoidal inputs, as defined in Eq. (5), can also be 
derived in terms of the integrator error parameters ej and k [3]. From this result the following 
asymptotic formulas are obtained for e\f, the fractional error in transfer function gain, and e\, 
the transfer function phase error: 


soo 



( 8 ) 


ill coXe k 

For k odd, e M = (-1) — j(coh) , e A 

co + A 

J m e i r t sk 

Fork even, e M = -(-1) “ 2 * e A 

co + A 


/t+l 2 

-r— flJ e, i. 

<- l > 2 • 
co + A 


(-D 


y ooAe 


2 ,2 

oo + A 


' (0)h) k 


coh « 1 


oo/t « 1 


(9) 


Here the errors are proportional to ej ( coh )*. Comparable asymptotic formulas for cm and ca 
can be derived for digital simulation of a second-order system with transfer function H(s) = 
1/(^2 + 2 t^cOnS + (On 2 ) [3]. Again, the gain and phase errors are proportional to <?/(« oh) k . 

The transfer function H(s) for any order linear system with both real and complex roots 
can be represented as the product of first and second-order transfer functions. In this case it can 
be shown that the asymptotic formulas for the digital system transfer function gain and phase 
eiTO rs is simply the sum of the individual first and second-order subsystem gain and phase 
errors, respectively, for predictor and predictor-corrector methods of the type shown in Table 1. 
If we simulate a flexible structure with a given integration method, this permits us to compute 
the linearized system gain and phase errors at the frequency (0 for any input-output pair as a 
function of integration step size h. In view of the reemerging popularity of frequency-domain 
methods for designing multiple input/multiple output control systems, this is a quite useful 
result. It permits us to estimate ahead of time for a given step size and integration method 
whether the simulation errors will be satisfactorily small. Conversely, for a given transfer 
function accuracy requirement, it allows us to compute the maximum allowable step size h for 
the simulation. 

It should be noted that the methodology outlined above for determining characteristic 
root and transfer function errors for any order of linearized system from the simple integrator 
model given by Eq. (6) does not work in the case of multiple-pass, single step methods such as 
Runge-Kutta. This is because the results of individual passes within a single step in such 
methods depend on the particular form of the system transfer function. Asymptotic formulas 
for the root error parameters e\, eo, and can, of course, be derived separately for RK-2, RK- 
3, RK-4, and variations of these methods [3]. 


3. Modified Euler Integration Algorithms 

In this section we describe some modifications of simple Euler integration which have 
potential advantages over conventional integration methods such as those listed in Table 1. First 
we introduce the concept of state variables defined at both integer and half integer sample times. 
Assume that the simulation of a mechanical degree of freedom with position state *, velocity 
state y, and acceleration a involves integrating the following simple state equations: 
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Next assume that successive data points are defined at integer time samples in representing the 
acceleration a and position x, and at half-integer sample times in representing the velocity y. 
The following modified Euler algorithms can then be used for integration: 


y n*in Vl* + h a n ’ X n* 1 X n + ^ n*m 01) 

The basic concept behind this modification of standard Euler integration is very simple; instead 
of the using the state variable derivative defined at the beginning of the integration step, the 
method uses a state variable derivative defined halfway through the step. For this algorithm it is 
easy to show that the integrator error coefficient defined in Eq. (6) is given by e f = 1/24 and the 
order of the method is k = 2. Thus the accuracy of this single-pass algorithm is twice that of the 
two-pass AM-2 predictor-corrector. However, the algorithm does require that the velocity 
states be defined at half-integer sample times. 

Let us apply this modified Euler method to the second-order system represented by Eq. 
(2) for the generalized coordinate x. We can replace Eq. (2) by the following two state 
equations: 

y = ® n i<t>(t)-x-2Cy] , x = 0 ) n y (12) 

By analogy with Eq. (11) the modified Euler difference equations become: 


n + l /2 


■ y 


n- 1/2 


+ " X n ~ 2 ^y ’n > ’ *„♦! = hy n 


+ 1/2 


(13) 


Since y„ is not explicitly computed, it is necessary to substitute an estimate y'„ in the damping 
term on the right side of the y n + 1/2 equation. There are many ways in which the y' n estimate can 
be computed. In Table 2 we list four candidate methods. 


Table 2. Methods for Estimating the Velocity y n in Modified Euler Integration 


Method 

Formula for the Estimate. 

1 . Average of y „ and y 

° -'n+l/2 J n- 1/2 

* _ y n+\/2 + y n-\/2 
y n 2 

2. Extrapolation using y n ^ and y n 

• 3 1 

y =^Ly -«-y 

J n 2 n -W 2 J n ‘ za 

3 . Integration using y n ,andy n 

- 7 . 3 

4. Estimate based on y 

J n-l/2 
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The estimate for y„ in the first method is simply based on averaging y n+l/2 and y n _ m . 
This is equivalent to utilizing trapezoidal integration for the damping term. Although this means 
that y n+ 1/2 now appears on both sides of the difference equation in (13), for the linear system 
considered here it is possible to solve explicitly fory„ + i/ 2 , as we will see in the next section. In 
the second method the estimate for y n is based on a linear extrapolation from y n . m and y„-3/2- 
This is equivalent to using trapezoidal integration for the damping term. Since y n+ 1/2 now 
appears only on the left side of the difference equation in (13), this method can be used in the 
simulation of equations where dy/dt is a nonlinear function of y. This is also true for the third 
and fourth methods. The third is based on a second-order predictor integration over the interval 
h/2, starting with y n -\H and using dy/dt at the n - 1 and n - 2 intervals. This is equivalent to 
estimating y n with quadratic extrapolation based on y n -U2, yn-Vl and y n -5n ■ In the fourth 
method we simply use y n .\n as our estimate for y„. This is equivalent to Euler integration for 
the damping term. 

4. Modified Euler Integration with Trapezoidal Damping 

We have seen in Table 2 that the velocity estimate y' n for the modified Euler difference 
equations in (13) can be based on the average of y n +in and y n -m- Thus 


y* m 


uz + y n . 


in 


(14) 


As noted earlier, this is equivalent to utilizing trapezoidal integration for the damping term. 
Although this means that y/ 1 + 1/2 now appears on both sides of the difference equation, for the 
linear system considered here it is possible to solve explicitly fory n+ i/ 2 . In this way we obtain 
the following equations: 


where 


>W = C 2^n' 




Jt , = + 

n +1 n 


(Ohy 


n-\fl 


c x 


l- Coy h coh 

* n Q _ n 

1 + Cco n h 2 1 + Ca) n h 


(15) 


(16) 


From the method of z transforms applied to Eqs. (15) and (16) we obtain the following 
asymptotic formulas for the frequency and damping ratio errors of the digital simulation [4]: 


e = 

CO 






l + 4£ 2 -8£ 4 2,2 

^~ 0 ) n h ’ 

24(1- C 2 ) 






1) 


co h « 1 

n 


The transfer function gain and phase errors are given approximately by 


(17) 


SO 3 


(18) 


Fractional 
gain error 


l//*l t 

Ml ' 1 ” eM ~ 


£t» 

-1 

■ 

! 1 


’ CO 
n 






: 1 

2 

r 3 

<N 

t 

+ 

M 


L®"1 


(CO h), 


coh « 1 


2CV 


Phase 


error 

- C A = 

2 

'-•t 

2 

2£cy 



+ 



CO 

n 


CO 

n 


(ah). 


coh « 1 


(19) 


The characteristic root errors in Eq. (17) and the transfer function gain and phase errors in Eqs. 
(18) and (19) are comparable with those for AM-2 integration for the same step size h [3]. Yet 
AM-2 is a two-pass method whereas the modified Euler with trapezoidal damping, as used here, 
is a single-pass method. Thus it will take only half as long to execute as AM-2 while producing 
comparable accuracy. Its accuracy is approximately 5 times better than the accuracy of AB-2 
integration when applied to the same second-order system. 

The accuracy of modified Euler integration when applied to a linear second-order system 
can be further improved by the technique of root matching, which was originally employed by 
Fowler to improve the performance of conventional Euler integration [5]. By taking the z 
transform of Eqs. (15) and (16) we can obtain exact analytic formulas for the undamped natural 
frequency (On* and damping ratio £* in terms of co n and £ From these formulas we can solve 
for (On and £ in terms of co„* and £*. If in these formulas we then replace (On and £ by co n ' and 
£> respectively, and (On* and £* by co n and £ respectively, we obtain the following [4]: 


co 



2 cos ( (Oji-J 1 - £ 2 ) 
cosh(£fU n /r) 


( 20 ) 


? = 


tanh (£n> A) 


( 21 ) 


If (On and £' from these formulas are used instead of co n and £ in Eqs. (15) and (16), the 
resulting digital simulation will exhibit (On* and £* values which exactly match the (O n and £of 
the continuous system being simulated. For a given step size h the co n \ £', C\ and C 2 can be 
precomputed, so that each integration step in simulating the second-order system only requires 3 
multiplies and 2 adds, as before. Now the charcteristic roots of the digital simulation will be 
exactly equal to those of the continuous system, regardless of the integration step size h. The 
approximate formulas for the transfer function gain and phase errors are given by [4]: 


5(94 



( 22 ) 


1 / LN 2 

e U B 24 {0)h) ' 


LLU 9 

(to hy , (oh« 1 

6o) n 


Note that the fractional error in gain, e M , is completely independent of the damping ratio C, and 
the phase error e A approaches zero as £ approaches zero. Thus our modified Euler algorithm 
with root matching will be especially effective in simulating lightly-damped second-order 
systems as will be the case in structural modes. This is illustrated in Figure 1, where gain and 
ph^rs Jfrequency for a second-order system with {=0.01 am plotted. Because of the 
sharp resonant peak in gain and the extremely rapid change in phase as 0) passes through a n , it 
is very critical that both the natural frequency and damping ratio of the digital simulation match 
that of the continuous system. The table at the bottom of the figure shows die transfer function 
errors for input frequencies in the vicinity of (O n for the specific case of (O n h = 0.5, which 
corresponds to only 2 integration steps per radian or 12.57 steps per cycle. Shown in the table 
are the gain and phase errors based on both an exact calculation from the system z transform, 
H*(ej(dh^ as W ell as the approximate formulas of Eq. (31). Note how closely the approximate 
caculations agree with the exact, even for the example here for which ah »0.5. 

Until now we have only analyzed the dynamic performance of the modified Euler method 
in the frequency domain. This has been accomplished by examining the gam and phase errors 
of the transfer function for sinusoidal inputs.We now consider the errors in computed response 
of the second-order system to a unit-step input. Figure 2 shows the errors which result when 
using RK-2 integration (Heun’s method); modified Euler with trapezoidal integration for the 
damping term, i.e., Eqs. (15) and (16); and modified Euler with root matching, i.e., On and C 
from Eqs. (20) and (21) substituted for On and £ in Eqs. (15) and (16). For the example m the 
figure the damping ratio £ = 0.707 and the integration step size is given by a n h = 0.5. The 
results show that the RK-2 errors are 4 to 10 times larger than the modified Euler errors. It 
should also be noted that RK-2 is a two-pass method, that is, it requires two evaluations of the 
state-variable derivatives per integration step. It follows that RK-2 will take approximately 
twice as long to execute per integration step as the single-pass modified Euler methods To 
provide the same output integration frame rate in real time the RK-2 method will therefore 
require twice the mathematical step size h in comparison with the modified Euler methods 
considered here. This will further increase by a factor of 4 the RK-2 errors relative to the 
modified Euler errors in Figure 2. 


The modified Euler results shown in Figure 2 were obtained using an initial step of hfl 
in integrating dydt to obtain y. After one integration step this provides the calculation of yi/2 
starting with the initial condition y 0 . The step size is taken as h for all subsequent dyldt 
integration steps. This results in successive velocity values representing y at half-integer step 
times, consistant with the concept introduced in the beginning of this section. 
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Fractional Gain Error 

Phase Error (radians) 


oVCDn 

Exact 

Eq. (31) 

Exact 

Eq. (31) 

co n h — 0.5 

0.7 

0.01040 

0.01021 

-0.000296 

-0.000292 


0.9 

0.01727 

0.01888 

-0.000381 

-0.000375 

12.57 steps 

1.0 

0.02137 

0.02083 

-0.000424 

-0.000317 

per cycle) 

1.1 

0.02592 

0.02521 

-0.000467 

-0.000458 


1.4 

0.04240 

0.04083 

-0.000595 

-0.000583 


Figure 1 . Frequency response of lightly-damped second-order system using modified Euler 
integration with root matching, C 0 nh = 0.5. 
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Figure 2. Unit step response errors in simulating a second-order system with damping ratio 
£=0.707, integration step size given by (Onh - 0.5. 


5. Performance of Other Versions of Modified Euler Integration 


In this section we present the asymptotic formulas for characteristic root and transfer 
function errors when modified Euler integration is used to simulate a second-order system with 
methods 2, 3, or 4 in Table 2 utilized to calculate the velocity estimate y' n in Eq. (13). For 
method 2, which is equivalent to AB-2 integration for the damping term, the following results 
are obtained for the fractional error in root frequency, and e;, the damping ratio error [6]: 


e 


CO 


1-32^ + 40 C 
24(1- £ 2 ) 



e e s 


11C-20C 3 , . ,2 

(co h) 

24 


(O h « 1 

n 


(23) 


These errors are significantly less than the errors when AB-2 is used for all integrations. For 
method 3 in Table 2, which uses a second-order predictor integration algorithm to compute y n , 
the following asymptotic formulas are obtained for the root frequency and damping errors. 


0-40, lx2 






co h « 1 

n 


(24) 


The transfer function gain and phase errors are given by 
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In both Eqs. (24) and (25) the errors are a factor of two smaller than the corresponding errors 
when AM-2 is used for all integrations. In addition, the AM-2 algorithm is a two-pass method 
which will therefore take twice as long to execute on a given computer. For method 4 in Table 
2, which is equivalent to using Euler integration for the damping term, the following formulas 
are obtained for the characteristic root and transfer function errors [4]: 



Note that the errors are all proportional to the first power of the step size h. This is because of 
the first-order Euler algorithm used for integration of the damping term. For £ = 0, however, 
the first-order errors in Eqs. 25) and (26) vanish, meaning that the errors become second-oider 
in h. This is to be expected, since the conventional Euler integration plays no role when £ = 0. 
In fact it can be shown that when £ = 0, the digital solution will have zero damping regardless 
of the step size h. & 

When method 2,3, or 4 in Table 1 (or any other explicit method) is used to provide the 
estimate y n for the velocity state, the modified Euler method can be used as the algorithm for 
integrating the nonlinear state equations represented by (1). The vector difference equations 
become the following: 


: + mq$Ciq n ,q) - K(qJ) , + A* 


We now turn to a consideration of integration algorithm stability. 
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6. Stability of Integration Methods 

It has already been pointed out that the stability of numerical integration algorithms 
becomes an important consideration when the flexible structure is modeled by discretization. 
This is because the discretized model will contain high frequency modes which are unimportant 
in the simulation but can cause numerical instabilities for reasonable integration step sizes. For 
a given integration method the stability boundary in the Xh plane can be obtained by considering 
a simulation of the linear system with transfer function H(s) = l/(s-A). From the difference 
equation the z transform, H*(z), is obtained. The stability boundary is defined by the Xh values 
for which the denominator of H*{z) vanishes when Izl = 1. These Xh values can be obtained by 
letting z=d e in the denominator of H*(z) and solving for Xh for 0 values ranging between 0 
and jc. When this is done for the AB predictor methods, the stability regions plotted in Figure 3 
are ob tain ed. The regions are symmetric with respect to the real axis so that only the upper half 
plane is shown. For any values of Xh lying outside the boundaries the digital simulation will be 
unstable. In the case of both AB-3 and AB-4 the boundary crosses over into the right half 
plane This means that a continuous system with roots on the imaginary axis which correspond 
to undamped transients can exhibit stable transients in the digital solution. Put another way, it 
means that AB-3 and AB-4 solutions will exhibit more damping than the continuous system 
being simulated. This is actually desirable in the case of the high frequency modes which are 
not of interest in a given simulation. On the other hand the AB predictor methods do not have 
particularly large stability regions and therefore do not permit very large integration step sizes h 
compared with the reciprocal magnitude, 1/1 Al, of the largest eigenvalues in the simulation. 



Figure 3. Stability boundaries for AB predictor integration. 


In Figure 4 the stability boundaries are shown for the the two-pass AM predictor-corrector 
methods. Although the boundaries are considerably larger than those for the AB methods, it 
must be remembered that the AM algorithms will take twice as long to execute. Thus the 
boundaries should be reduced by a factor of two for a valid comparison with AB-2. When this 
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is done, the AM-2 and 3 boundaries actually fall inside the AB-2 and 3 boundaries, although the 
AM-4 boundary still lies outside the AB-4 boundary. In all cases the higher-order algorithms 
exhibit less stability and are therefore unlikely to be candidates for simulating flexible structures 



Figure 4. Stability boundaries for two-pass AM predictor-corrector integration. 

For comparison purposes the stability boundaries for RK-2, 3 and 4 are shown in Figure 5. 
We recall that these algorithms require 2, 3 and 4 passes, respectively, through the state 
equations per integration step. Thus for proper comparison with single-pass methods the 
boundaries shown should be reduced by factors of 2, 3 and 4, respectively. When this is done, 
the RK-2 boundary roughly matches the AB-2 boundary, while the RK-3 and RK-4 boundaries 
still fall outside the AB-3 and 4 boundaries, respectively. 



Figure 5. Stability boundaries for Runge-Kutta integration methods. 
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F inall y, in Figure 6 are shown the stability boundaries for various modified Euler 
methods, as described in Sections 4 and 5. The trapezoidal damping case corresponds to 
method 1 in Table 2, the Euler damping case to method 4, the AB-2 damping case to method 2, 
and the predictor damping case to method 3. Also shown for comparison purposes in Figure 6 
are the stability boundaries for AB-2, AM-2 and RK-2, as presented earlier in Figures 3, 4 and 
5, respectively. The AM-2 and RK-2 stability boundaries have been reduced by a factor of two 
to reflect the two passes per integration step required in the implementation of these methods. 
Note that all four of the Modified Euler methods in Figure 6 have stability regions which permit 
values of \Xh\ up to 2 for lightly damped transients, e.g., eigenvalues near the imaginary axis. 
In this regard the methods are considerably superior to the AB-2, AM-2 and RK-2 algorithms 
and should perform especially well in the simulation of flexible structures. 

It should also be noted that the modified Euler methods are ideally suited for real-time 
simulation in that they do not require inputs prior to their occurence in real time. For example, 
if p(t) in Eq. (1) is a real time input, then the single-pass modified Euler algorithm of Eq. (28) 
only requires F n at the beginning of the nth integration step. On the other hand, the AM 
predictor-corrector algorithms require F„+i at the start of the second pass for the nth integration 
step, and F n+ i is not yet available in real time. There is, however, a modified version of the 
AM-2 predictor method which is compatible with real-time inputs [6]. The AB predictor 
methods are also compatible with real time inputs, and there are versions of RK-2 and RK-3 
which permit real-time inputs [3]. RK-4 is not compatible with real-time inputs, since it 
requires F n +in at the beginning of the second pass and F n+ i at the start of the fourth pass, in 
both cases prior to their availability in real time. 



Figure 6. Stability boundaries for modified Euler integration methods. 

ill 








7. Conclusions 


In this paper we have considered the dynamic performance of integration methods in the 
context of simulating flexible structures. In terms of both characteristic root errors and transfer 
function errors, both important in such simulations, we have compared the performance of 
traditional integration methods with various versions of modified Euler integration . We have 
shown that modified Euler integration is especially effective in simulating lightly-damped 
structural modes. We have also shown that the modified Euler methods have very favorable 
stabilty boundaries in the Xh plane with respect to requirements in the simulation of lightly- 
damped modes. This is especially significant when a flexible structure is modeled by 
discretization as opposed to normal coordinates, since it will allow larger integration step sizes 
before the solution goes unstable due to the presence of higher modes which are unimportant to 
the simulation. 


References 

1 . Gilbert, E.G., "Dynamic Error Analysis of Digital and Combined Digital- Analog 
Systems," Simulation, vol. 6, no. 4, April 1966, pp 241-257. 

2. Benyon, P.K., "A Review of Numerical Methods for Digital Simulation," Simulation vol. 
11, no. 4, Nov. 1968. 


3 . Howe, R.M., "Transfer Function and Characteristic Root Errors for Fixed-Step Integration 
Algorithms," Transactions of the Society for Computer Simulation, vol. 2, no. 4 Dec 
1985, pp 293-320. 


4. Howe, R.M., "Simulation of Linear Systems Using Modified Euler Integration Methods," 
to appear in Transactions of the Society for Computer Simulation. 

5. Fowler, M.E., "A New Numerical Method for Simulation," Simulation, vol. 4, no. 5 Mav 

1965, pp 324-330. ’ * 

6. Howe, R.M., "The Role of Modified Euler Integration in Real-Time Simulation," 
Proceedings of the Conference on Aerospace Simulation II, San Diego, 1986; pp 263-275 
The Society for Computer Simulation, P.O. Box 17900, San Diego, CA 921 17. 




THIS PAGE INTENTIONALLY LEFT BLANK 


OPTICAL PROCESSING FOR DISTRIBUTED SENSORS IN 
CONTROL OF FLEXIBLE SPACECRAFT 


- Id 
>- 06 
06 O 
U Z 
Z M 
O -I 
O -J 

H « 
z o 
o 


o z 
z 

O Q 

z z 
> < 


06 

Id 

I- I 
Z O in 
Id Z CM 
O < CM 

k in 
z ca i 
o in 

06 mj VO 
< OVD 
Ld 06 CO 
(A h CM 
Id Z 
O* O 


06 

* 

>- 

>■ 

Id 

l- 


* 

-1 

b. - 

Z 

(/> 

a 

< z 

CJ 

a 

z 

a o 

-1 < 
id co 

.3 

O K- 

ui a. 
u z 


in 

-i 
Z Id 

2 < 
« z 

< o 

Z H 

in z 


cn o. 
< in 


id 

H 

a 

o 

u. 


o 

H 

h* 

£ 


id 

Id 


514 


PRECIDrMG Pa nr r>, uu 

RU.fOK MOT FJUWEi 


Workshop on Computational Aspects in the 
Control of Flexible Systems 
Williamsburg, VA 



